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ABSTRACT 



An examination and a comparison of the relative merits of 
the finite element and Fourier series methods of solving radially 
loaded circular ring problems are made. The procedure employed 
to evaluate the two methods is to use each method to solve for three 
different load conditions and to compare the performance of the two 
methods on the basis of accuracy, ease of usage, and equipment 
required. The results indicate a satisfactory accuracy for both 
methods under most conditions. The Fourier series method is 
superior for solving problems with a distributed load condition. 

The finite element method is superior for solving problems with 



concentrated loads. 
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CHAPTER I 



INTRODUCTION 



The structural analysis of many aerospace vehicles, such as 
rockets, re -entry bodies , aircraft, space capsules, etc., is 
frequently accomplished using either the Fourier series method 
or the finite element method in conjunction with a digital com- 
puter. A considerable amount of effort has been devoted to the 
development and refinement of these two methods, but little 
attention has been given to a direct comparison of the methods. 
Fundamental questions regarding the similarities, differences, 
relative accuracy, ease of usage, and advantages or disadvantages 
of each of the methods have not been answered as yet. The ob- 
jective of this study was to seek answers to some of these funda- 
mental questions. 

The procedure employed to accomplish the objective was to 
apply the two methods of analysis to a portion of a typical aero- 
space structure and to investigate the relative merits of the two 
methods. The portion selected for analysis was a thin circular 
ring with a rectangular cross section. The ring was selected 
because it represents, in a simplified manner, the circular 
cylinders that are frequently found in aerospace vehicles. Three 
different load conditions were applied to the ring for each of the 
two methods of analysis so that the effects of the load distri- 
butions on the performance of the two methods could be investi- 
gated. The extent to which the analyses were carried out in 
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terms of numbers of significant figures carried, numbers of 
terms evaluated in the Fourier series method, number of 
elements used, and use of double precision techniques in the 
FORTRAN computer programs, was greater than that normally 
used in common engineering practice. The reason for such 
depth was to bring out any subtleties peculiar to either of the 
two methods that might not appear if normal engineering accu- 
racy were used. 

The author gratefully acknowledges the guidance and 
assistance of Professor Robert E. Ball in the preparation 
of this thesis. 
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CHAPTER II 



DERIVATION OF THE GOVERNING EQUATIONS 
FOR THE FOURIER SER IES A NA L YSIS 

Consider the thin ring shown in Figure 1 where a is the radius 
of the ring middle surface; b, the ring width; and t, the thickness. 
The angle ([) is the generator of the ring and is positive in a counter 
clockwise direction. All points in the ring are located by the 
coordinates 0 and z, where z originates at the middle surface and 
is positive in the outward radial direction. The normal force N^j 
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the shear force Q^, the bending moment M^, and the applied 
radial loading P r are shown in the positive direction. The dis- 
placements, also shown in the positive direction in Figure 1, are 
w in the radial direction and v in the tangential direction. 

Assuming that the ring retains its circular shape after de- 
formation, the equations of equilibrium can be written as 

aQ<K = 0 ( 2-la ) 

a 4> 



and 




(2-lb) 



(2-lc) 



Using Equation (2-la) to eliminate from Equations (2-lb) and 
(2-lc) leads to the two equilibrium equations 



^ d M 4> 

A 4> <A 4) (2-2a) 



aN^- a 2 P r + 



fa* . 



= o 



♦ " -r a _ (2-2b) 

The circumferential strain, £ ^ , at any distance z from 



the middle surface, is given in Reference 1 as 

1 2 






i 



V 



z. i a tu u) 

+ 



aAcji <\ (a. + 2 .) a+z. 



(2-3) 



This relationship is based upon the assumption that points on a 
normal to the middle surface prior to deformation remain on the 
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normal after deformation and that the thickness of the ring 
remains constant. Substituting Equation (2-3) into the plane 
stress form of the elastic constitutive law leads to 



<r, = E 




z i IQ 

a. (a + z) a + z 



(2-4) 



where (f ^ is the circumferential stress at the location (<|), z) 
and E is the modulus of elasticity. 

The normal force and bending moment are given by 



t 







and 




Z 



Substituting Equation (2-4) into Equations (2-5a)and (2-5b), 
integrating and applying the limits yield 



( 2-5b) 



N* = 



_D_ 

a 







K 


d Z u> ] 


t- uu 




a 3 


+ UL) 



and 



K . _ K_ 

a*- 









d fy 2 - 

where D - Et and KzEt^/12 . * 



(2 -6a) 



(2- 6b) 



*In order to use this form of D and K the natural logarithmic 
terms that result from the integration of Equations (2-5a) and 
( 2 - 5b ) must be expanded in a series in terms of t/2a where terms 
to the fifth power and above are dropped. 
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The problem has now been reduced to one with four unknowns, 
v, w, N^, and M^, and four equations; two of the equations are 
obtained from the conditions of equilibrium 



dM* 
a — - 
d<t> 



d4> 



= o 



(2- 2a) 






A <t> x 



(2-2b) 



and two are obtained from the application of the elastic consti- 
tutive law and the strain-displacement relationship 



N*= 



D 

a 2 - 



ci\T 

d<> 



4- UU 



K 



d (p 2- 



4- UU 



( 2-6a) 






K 

a}- 



d^uu 

TT 2 - 



+ 0 ) 



( 2 - 6b ) 



Taking the appropriate derivatives with respect to (j) of Equations 
(2-6a) and (2-6b) and substituting the results into Equations (Z-2a) 
and (2-2b) yield 



d T d V 

d<J) L d<t> 



4- 10 



o 



(2-7a) 



and 



dvr 



+- 0 ) + 




-t-Z 



d Z tu 
<1 (p 1 



LL) 




D 



( 2- 7b) 



O O t 

where k = t /12a . Equations (2-7a) and (2-7b) are the governing 
differential equations for a thin ring under the assumptions that 
have been made. 
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CHAPTER IH 



DEVELOPMENT OF THE FOURIER SERIES SOLUTION 
TO THE GOVERNING EQUATIONS 

The starting point in the development of the Fourier series 
solution to Equations (2-7a) and (2-7b) is to assume that the 
applied load, P r , is a general function of 0, symmetrically dis- 
tributed about two perpendicular planes, one through (j) = 0 an d 
0 =: tT, and one through 0=^1T/2. Thus, the load can be repre- 
sented by the Fourier cosine series * 



for the ring under the load distribution considered, is also written 
in a Fourier cosine series as 



w> - o, z, *+ 



* The symmetry conditions on P r lead to the elimination 
of all odd values of m, i. e. , P m = 0 for all odd values of m. 




(3-1) 



TK - 0 , 2,4 

where the coefficients P m are given by 




o 



( 3-2a) 



and 




O 



( 3-2b) 



The radial displacement, w, which is an even function of (j) 



oo 




(3-3) 
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The tangential displacement v is an odd function of (J) and is 
given by 



oo 

V = £ Mn m 4> 0-4) 

The governing differential equations are used to obtain w^ 

and v as a function of P„. For the axisymmetric mode, 
m m J 

where m = 0, the substitution of Equations (3-1) and (3-3) into 
Equation (2-7b) yields 



uu 



o 



a 1 P 0 
0 (1 +Sk) 



(3-5) 



For the general case of m >0, substitution of Equations (3-3) 
and (3-4) into Equation (2-7a) results in 



\r m + m - o 



(3-6) 



Substituting Equations (3-1), (3-3), and (3-6) into Equation (2-7b) 
gives 

l 

Jk D ( YA l --i) Z 

Therefore, according to Equation (3-6) 

m 

( 

Substituting these expressions for w m and into Equations (3-3) 
and (3-4) leads to. 



\Tw\ = — 



r 



A D 



m 




a 1 P 0 a 1 o? P m cos w\ 4? 

D(HA') + AD 1 

tv\- 
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( 3 - 7b ) 



V = 



a 



Jk D 



CO 

z 



p 



m 



w(w x -iY 



m - 2..^ 

All other quantities of interest are determined by using these 
expressions for w and v in the appropriate equations. For 
example, according to Equations (2-6a) and (2-6b), the normal 
force is given by 



N*=aP 0 - a I 



W\ 



COS W» ($> 



*<$> 

m = 

And the bending moment is 






= - 






CO 



1-mt 



- & L 



Pyy^ COS W\ (p 



vv\ = 



W\^“ - 1 



(3-8) 



(3-9) 



To permit a direct comparison between the Fourier 
series and finite element analyses of the ring, w, N^, and 
are selected for evaluation by both methods at the points (() ~ 0, 
30, 60, and 90 degrees. In Chapter IV, P^, and the specific 
equations for w, N^, and are developed and evaluated at 
these values of for three different load conditions. 
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CHAPTER IV 



APPLICATION AND RESULTS OF THE FOURIER SERIES SOLUTION 

The three loads selected for use in this study are shown in 
Figure 2. All three loads are radial loads, directed inward, and 
are symmetrically distributed about two planes, one through 0=0 
and 0 = IT and one through 0 = "tlT/2. Case 1 is a "line" or "knife 
edge" load of P pounds distributed over the width of the ring at 0 — 0 
and 0 = TT* 




Case 


i 


OC z 


0 degrees 


Case 


2 


<K T 


30 degrees 


Case 


3 


K = 


60 degrees 



FIGURE 2 

LOADS CONSIDERED 
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Case 2 and Case 3 loads are uniform pressure loads, each with a 
total value of P pounds covering the areas shown in Figure 2* The 
uniform loads are initially represented as P pounds distributed 
over the angle 2 e < so that the surface load P r is given by 

P 

— (4-1) 

2. e< A u 

Thus, substituting Equation (4-1) into Equation (3-1) gives 



_P 

a b 



-z C03 4> 

m = o, 



Accordingly Equation (3-2a) becomes 



Po 



so that 



P 

2®<a b 





P 

tTab 



(4-2) 



Siroiliarly, substituting Equation (4-1) into Equation (3-2b), inte- 
grating and applying the limits yield 



• 1 

P =- - s \ 'A rA (TT— °<)| 

m IT a b vY\<=< 

Therefore 

j. zP sm m=< 

M IT a b vy\ cx. 

for m even. For odd values of m 



(4-3) 




= o 



as originally noted in Equation (3-1). 
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Applying Equation (4-2) and (4-3) to Equations (3-7a), (3-8) 



and (3-9) leads to the Fourier series solution for w, 



N and M. 
<P <1> 



UU - - 



C.P 



2 a P 






® S w\ot c os w\ (J> 

^ YY\ (m 1 - !) 2 - 

m = z, *t 



( 4 - 4a ) 



N 



4> 



P ZP 2 ? cos m 4 > 

+ 7 

irV> trb w\ ( 

Wi = 



( 4 - 4b ) 



and 



u _ _ . 2jd 

ffkfl+A) 1Tb 

The non-dimensional forms of these three equations, denoted by 
w, and M^, and obtained by using A = t^/ 12a^, D = Et, and 
I = bt^/ 12, are 



2? $> W\ VY1<=< C O S> VY> (t> 

£ 1 ( 4-4c ) 

YYifm 1 - lV 

m = z, h 



uoEl 

UL) - 

a>P 



I 

IT a 1 A(l + i) 



z. PS sm m<=< c os w\(t> 

IT ~ °< m ( yv\ 2 - iV ~ 

M 



(4-5a) 



N <f) = 



kN* 



1 ® S\v\w\oc COS vn ((> 

+ L 



ft *“ =< W\ ( W1 *• -l) 



(4-5b) 



and 



_ Vs Wl Jj I L f StV\YY\c< COS d, 

M 3 £ = q S — 

<\P ITa*- A(l+i^ IT ^w\(yaI-±) 

m = z,*f 



(4-5c) 
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These expressions are applicable only to Case 2 and 3 loads. 
For Case 1 load, c< — 0. Thus, taking the limit of Equations 
(4-5a), (4-5b) and (4-5c) as c* — >. 0 leads to the solution for 
Case 1 load 



x 



% 

IT 



z - 

m - 



cos m 4> 
(m z - 1) x 



(4-6a) 



N $ = - 



and 



IT 



-) 00 
l 



co^ m d) 



ir va z -1 

m -x,H- 



(4-6b) 



M*=- 



7 oo 



COS YV\ d) 



lTa z Ml+Jk) TT YYl z -t 

m = i, h 



(4-6c ) 



The term - I/U A (1+ k) that appears in w and will now 
be examined for its contribution and significance. The expressions 

dU JiV^ and j U - MIaM 

a iAE b 2. E I 



give the axial strain energy and the bending strain energy respec- 
tively of an element V>ta(d(j)) of the ring. They show that the axial 
strain energy is proportional to 1/A and the bending strain energy 
is proportional to 1/1. Using 1/A and 1/1 leads to 

U a /U^^ 1/A. In Reference 2, the axial strain energy for such a 
curved element is shown to be very small compared to the 
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beading strain energy of the element provided that t/a is small. * 

Thus, the term - I/1T A (1 + k) in w and M, is related to the 

extension of the middle surface of the ring and should be small. 

If the middle surface is considered inextensional, the axial 

strain energy of an element of the ring is zero, and the term 
o 

-I/1T a 4 * A(l-t-k), which is proportional to 1^/11^, does not 
appear in the expressions for w and M^. In order to determine 
the magnitude of the contribution of the axial strain energy to the 
solution, a t/a ratio must be selected. As an example, a ratio 
of t/a = 1/10 is chosen and w and are calculated with and 
without the term containing 1/A at (J) = 0, 30, 60 and 90 degrees 
for all three load conditions. The results of these calculations 
are found in Table I, where the values referred to as "EXT" 
include the -I/TT A (1 + k) term and the "NON-EXT" values do 
not. A s is evident from Table I the differences between the "EXT" 
and n NON-EXT n values are very small, or nonexistent for the 
accuracy retained, when t/a = 10. 

An important consideration in the application of the Fourier 
series method is the number of terms of the series required 
to give a desired accuracy in the solution, i. e. , how fast does 



*This statement assumes there are no distributed loads on 
the element, i. e. , only concentrated loads are considered. This 
thesis considers distributed loads, and thus the statement is not 
completely applicable. 
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TABLE I 



FOURIER SERIES EXT ENSIONA L AND NON - EXTENSIONAL RESULTS 





CASE 1 


LOAD 


CASE 2 LOAD 


CASE 3 


LOAD 






Ext. 


Non-ext. 


Ext. 


Non-ext. 


Ext. 


Non-ext, 




w 


-0. 0747 


-0.0744 


-0. 0599 


-0. 0596 


-0. 0289 


-0. 0287 


0° 




0 . 0 


0 . 0 


-0. 128 


-0. 128 


-0. 239 


-.0. 239 






0. 318 


0. 318 


0. 190 


0. 190 


0. 079 


0. 080 




w 


-0. 0337 


-0. 0334 


-0. 0289 


-0. 0287 


-0. 0152 


-0. 0149 


30° 




-0. 250 


-0. 250 


-0. 239 


-0. 239 


-0. 271 


-0. 271 




M 4 


0. 068 


0. 068 


0. 079 


0. 080 


0. 047 


0. 048 




w 


0.0361 


0. 0364 

I 


0. 0295 


0. 0298 


0. 0141 


0. 0143 


60° 




-0. 433 


-0. 433 


-0. 413 


-0. 413 


I -0. 358 


-0. 358 






-0. 115 


-0. 115 


-0. 095 


-0. 095 


-0. 040 


| -0. 040 




w 


0. 0680 


0. 0683 


0. 0571 


0. 0574 


0. 0295 


0. 0298 


90° 


N. 

<t> 


-0. 500 


-0. 500 


-0. 477 


-0. 477 


-0.413 


-0. 413 




M. 


-0. 182 


-0. 182 


-0. 159 


-0. 159 


-0. 095 


-0. 095 
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the series converge? This information for the inextensionai 
form of the solution is shown in Tables II, III and IV for the first 
4 terms of w, N^and at (J) r 0, 30, 60 and 90 degrees for all 
three load conditions. The number of terms required for con- 
vergence of a series to four decimal place accuracy for the 
displacement and three decimal place accuracy for the forces 
and bending moments is shown in the last column with the con- 
verged value of the series. 
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TABLE II 



FOURIER SERIES NON -EXT ENSIONA L CONVERGENCE 
FOR CASE 2 LOAD _ 



Term 


i 


2 


3 


4 ! Converged 




_JL 


m 


2 


4 


6 


8 


Value 


l 




2?m 


-0. 0707 


-0. 07.36 


-0. 0741 


-0. 0742 


-0. 0744 at m = 


14 ! 


0° 




-0. 106 


-0. 064 


-0. 045 


-0. 035 


0. 000 at id = 


534 






0. 212 


0. 255 


0. 273 


0. 283 


0. 318 at id = 


394 




L -' m 


-0. 0354 


-0. 0340 


-0; 0334 




-0. 0334 at id =r 


6 


3 0° 


3l_ 


-0. 212 


-0. 233 


-0. 252 


-0. 257 


-0.250 at id — 


34 




“ m 


0. 106 


0. 085 


0. 067 


0. 062 


0. 068 at id = 


52 




BVn 


0. 0354 


0. 0368 


0. 0363 


0. 0363 


0. 0364 at m = 


14 


60° 




-0. 424 


-0. 446 


-0. 427 


-0. 432 


- 0. 433 at id = 


24 






-0. 106 


-0. 127 


-0. 109 


-0. 114 


- 0. 1 1 5 at id = 


38 




§m 


0. 0707 


0. 0679 


0. 0684 


0. 0683 


0. 0683 at m : 


8 


o 

o 


S 

^ m 


-0. 531 


-0. 488 


-0. 506 


-0. 496 


-0.500 atm = 


26 




^ ID 


-0. 212 


-0. 170 


-0. 188 


-0. 178 


-0. 182 atm = 


38 



% 
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TABLE III 



FOURIER SERIES NON -EXTENSIONA L CONVERGENCE 

FOR CASE 2 LOAD 



Term 


1 


2 


3 


4 


Converged 


<t> 


m 


2 


4 


6 


8 


Value 




IXn 


-0. 058 


-0. 0547 


-0. 0597 


-0. 0596 


- 0. 0596 at m = 8 


0° 




-0. 143 


-0. 125 


-0. 125 


-0. 127 


-0.128 at m = 10 




Sift, 


0.175 


0. 193 


0. 193 


0. 191 


0.19 0 at m = 20 






-0. 0292 


-0. 0287 






-0. 0287 at no = 4 


30° 




-0. 231 


-0. 239 


-0. 239 


-0. 238 


-0.239 at m = 10 

• . i 






0. 088 


0. 079 


0. 079 


0. 080 


0. 080 at m = 14 






0. 0292 


0. 0298 






0. 0298 at m = 4 


o 

o 

o 




-0. 406 


-0. 415 


-0. 415 


-0. 414 


-0.413 at m = 44 






-0. 088 


-0. 097 


-0. 097 


-0.095 


-0.095 at m = 8 






0. 0585 


0. 0573 


0. 0573 


0. 0574 


0. 0574 at m = 8 


o 

o 

o 




-0. 494 


-0. 476 


-0. 476 


-0. 478 


- 0. 477 at m = 10 




Em™ 


-0. 175 


-0. 158 


-0. 158 


-0. 160 


- 0. 159 at m = 10 
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TABLE IV 



FOURIER SERIES NON -EXT ENSIONA L CONVERGENCE 

FOR CASE 3 LOAD 















) 














Converged 


Term 


1 


2 


3 


4 


Value | 

1 


° , 




2 


4 


6 


8 


t 

i 

» 

i 














j 


i 


' w 
- m 


-0. 0292 


-0. 0287 






-0. 0287 at m = 4: 


0 ° 




-0. 231 


-0. 239 


-0. 239 


-0. 238 


f 

- 0. 239 at m n 10! 






i 

0, 088 


0. 079 

j . 


0. 079 

1 


0. 080 


0. 080 at m = 14! 




w 

m 


-0. 0146 


-0. 0149 




t ' ' '■" l 

i 


“0. 0149 at m= 4 


30° 




-0. 274 


-0. 270 


-0. 270 


-0. 271 


-0.271 at mr 8 


l 


Lm 

m 


0. 044 


0. 048 


0. 048 


0. 048 


0. 048 at m = 14 




I™m 


0. 0146 


0. 0143 






0. 0143 atm= 4 


O 

o 

vO 


*-N m 


-0. 362 


-0. 358 






- 0. 358 at m = 4 


I 




-0. 044 


-0. 039 


-0. 039 


-0.040 


-0. 040 at m = 8- 




Z*rr, 


0. 0292 


0. 0298 






0. 0298 at m - 4 


o 

o 


T$ m 


-0. 406 


-0. 415 


-0. 415 


-0. 414 


-0. 413 at id = 44 


1 


Z-Mm 


-0. 088 


-0. 097 

1 


-0. 097 


-0.095 


-0. 095 at ra - 8 
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CHAPTER V 



DEVELOPMENT OF THE FINITE ELEMENT SOLUTION 

The finite element method is based upon the concept of 
replacing the actual continuous structure with a number of struc- 
tural elements of finite size, leading to an assembled structure. 
The requirements of equilibrium and continuity at the element 
joints or nodes lead to a set of simultaneous algebraic equations. 
There are two general methods for formulating these equations, 
the displacement method and the force method. The displace- 
ment method, also known as the stiffness or direct stiffness 
method, is the one used in this thesis to solve the ring problem. 
This method is discussed in detail in Reference 3. A brief ex- 
planation of the method is given here for continuity. 

In the displacement method the displacements and rotations 
at the nodes of the elements comprising the structure are taken 
as the unknown quantities. The matrix equation which relates 
the external forces to the displacements of an assembled structure 
is 



that act upon the nodes of the element. In a general three dimen- 
sional structure these may include bending and twisting moments 



symmetric matrix for the entire structure. This matrix is 
known as the stiffness matrix and is assembled from the individual 
stiffness matrices of the elements forming the structure. The 




(5-1) 



Matrix \ is the column matrix of nodal forces, i. e. , forces 
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entries in the element stiffness matrices are called stiffness 



influence coefficients. They are determined by displacing a 
node of the element by one unit in the direction of one. of. the 
degrees of freedom permitted for that node and calculating 
the nodal forces required to maintain that unit displacement and 
to prevent movement in any other degree of freedom for the 
nodes on the element. The amount of force required at a node 
to maintain this displacement state is one stiffness influence 
coefficient. The complete set of stiffness influence coefficients 
is determined by repeating this procedure for each degree of 
freedom at every node. The column matrix \u\ is the matrix of 
displacements and rotations which may take place at a node. The 
column matrix [x°j is a matrix of initial forces. They are deter- 
mined by clamping all of the nodes of the assembled structure 
and calculating the external force required at each node to 
maintain this zero displacement state under the applied load con- 
dition. Thus, loads distributed between nodes and concentrated 
loads not located at nodes can be treated. 

There are various techniques available for the manipulation 
and solution of Equation (5-1). The one used in this study is to 
transpose |x°] to the left-hand side of Equation (5-1). This leads 



the matrix equation is reduced in order by a consideration of the 
boundary conditions on the structure. If the boundary conditions 
specify zero displacement at some nodes, then those nodal dis- 



to a single column matrix of forces 




. Next, 



placements are no longer unknowns and are removed from 
along with the row of the equation and the column of 
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correspond to that zero displacement. The remaining matrix 
equation, which is referred to as the reduced matrix equation and 
is subscripted with /° , is 




The unknowns in Equation (5-2) are contained in the matrix of 
displacements ^u^J . Equation (5-2) is a set of simultaneous 
linear algebraic equations that is solved on a digital computer 
using the FORTRAN IV subroutine M DSIMQ, M a description of 
which appears in Appendix C. The outputs of M DSIMQ n are the 
nodal displacements i^Up| . 

In order to determine the nodal forces {x\ , the stiffness 
matrix equation is solved for each element. These equations, 
subscripted Pi , are 



i x 4 = 



(5-3) 



The unknowns in this equation are in j X^| since the displacements 
u ^ are obtained from the appropriate locations in \ u f\ anc * (uj. 
The matrix jXy|| is obtained by adding the matrix jx’jj] , which 
is determined by selecting values from appropriate locations in 
^ X°| , to the matrix jK^j |u^| . With iuy^ and jx^ known, 

all of the nodal displacements and forces have been determined, 
and the problem is solved since the forces and displacements at 
any interior location of an element can be determined from the 
applied load and nodal forces and displacements. 
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Applying the Unite element method to the ring problem under 
consideration requires the selection of the element geometry. 
Either straight or curved elements can be used. Curved elements 
are selected since they more naturally match the shape of the ring 
and result in displacements and forces that are readily interpreted 
without requiring a change in coordinates. The curved element 
used is drawn in Figure 3 with the nodal forces N , Q , N , Q , 

i 1 b b 

bending moments M^, displacements w^, v^, w^, v^ and 

rotations 0^, 0^ shown in the positive sense. The element covers 
the angle ft . A stiffness matrix for this element has been 
derived in Reference 3* based on the assumption that the thick- 
ness of the element, t, is small compared to the radius of 



❖A more detailed derivation and a comprehensive study of 
the curved element appear in Reference 4. 











FIGURE 3 



FINITE ELEMENT CONFIGURATION 
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curvature, a, so that only the strain energy of bending need be 
considered during the derivation. This assumption neglects the 
shear distortion energy, the energy resulting from any coupling 
between bending moment and normal force, the energy of axial 
deformation, and the displacements that correspond to these 
three forms of strain energy. The stiffness matrix equation from 
Reference 3 for the curved element can be given in a form corre- 
sponding to Equation (5-2) as 



N. ^ 


(v > 




N.° 1 


1 




1 




1 


Q 1 




- W 1 




V 


M i /a 




a0 l 




Mj°/a 


] N 2 


> = 7T M • 




■ + < 


N ° 
2 










o 


Q 2 




-w 2 






^4 ^ / 3 , 

\ ^ J 




* ae 2 




, M 2°/^ 



(5-4) 



where j^K^Jis the six by six matrix of stiffness influence 
coefficients for the curved element. The matrix K^j , each 
term of which is a function only of ft , will not be written out in 
functional form because it is quite lengthy. For this thesis a 
double precision digital computer subroutine, "RELM, 11 was 
written which computes ^K^jfor an input of ft . A copy of 
"RELM" is included in Appendix C. The stiffness matrix j^K^J 
computed by "RELM n for /3 = 7. 5, 15, and 30 deg rees is 

presented in Tables V, VI and VII respectively. 

With determined, can be assembled. Equation (5-1) 

is completed for the ring by using the loading conditions to obtain 
^X°| and the appropriate portions of |x( . Because of the 
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TABLE Y 





K) 


for 


/3 - 7. 5 degrees 




18, 665, 685 


1, 223, 063 


26, 675 


-18, 665, 640 


1, 223, 762 


-26, 721 


1, 223, 063 


85, 518 


2, 099 


-1, 223, 762 


74, 855 


-1, 400 


26, 675 


2, 099 


69 


26, 721 


1, 400 


-23 


-18, 665, 640 


-1, 223, 762 


-26, 721 


18, 665, 685 


-1, 223, 063 


26, 675 


1, 223, 762 


74, 855 


1, 400 


-1, 223, 063 


85, 518 


-2, 099 


-26, 721 


-1, 400 


-23 


26, 675 


-2, 099 


69 



TABLE VI 

K^] for 0=15 deg rees 



576, 897 


75, 862 


3, 306 


-576, 874 


76, 035 


. -3, 329 


75, 862 


10, 658 


524 


-76, 035 


9, 339 


-350 


3, 306 


524 


34 


-3, 329 


340 


-11 


-576, 874 


-76, 035 


-3, 329 


576, 897 


-75, 862 


3, 306 


76, 035 


9, 339 


350 


-75, 862 . 


10, 658 


-524 


-3, 329 


-350 


-11 


3, 306 


-524 


34 
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TABLE VII 



M for (% =30 degrees 



17, 243 


4, 598 


399 


-17, 232 


4, 640 


-411' 


4, 598 


1, 317 


130 


-4, 640 


1, 158 


-87 


399 


130 


17 


-411 


87 


-6 


-17, 323 


-4, 640 


-411 


17, 243 


-4, 598 


399 


4, 640 


1, 158 


87 


-4, 598 


1,317 


-130 


-411 


-87 


-6 


399 


-130 


17 
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symmetrical nature of the applied loads only one quadrant of the 
ring needs to be considered. The quadrant from (j) = 0 to 90 degrees 
is selected so that a direct comparison with the Fourier series 
solution for w, and can be made for ())= 0, 30, 60 and 

90 degrees . 



/ 
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CHAPTER VI 

APPLICATION AND RESULTS OF THE FINITE ELEMENT SOLUTION 



In order to apply the finite element solution to the ring 
quadrant (j) = 0 to 90 degrees, Equation (5-1) is assembled for the 
entire quadrant using Equation (5-4) for each of the elements com- 
prising the quadrant. ^KqJ matrices for ft — 7. 5, 15 and 30 de- 
grees are used in this thesis so that the quadrant is divided into 
12, 6 and 3 elements respectively for the solution of each of the 
three load conditions. This assembled equation is reduced and 
solved in the manner described in Chapter V. The process of 
assembling, reducing, and solving Equation (5-4) is best described 
by using an example. Case 2 load with ft = 30 degrees is selected 
as the example. 

Figure 4 shows the assembled quadrant and the node numbering 
sequence for the three elements in the quadrant. For three 




FIGURE 4 

3 ELEMENT QUADRANT FOR CASE 2 LOAD 
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elements ! k| is a 12 x 12 square symmetric matrix which is 
assembled using three of the 6x6 matrices from Table VII. 

When assembled for this example, Equation (5-1) appears as 



N, 



Qj = 0 

Mj/a 



n 2 = 0 
q 2 = 0 



) M 2 /a = o El k| 12 X 
\ ! a ^ ' 



12 



n 3 = 0 
q 3 = 0 

M 3 /a 

N„ 



Q 4 = 0 

M /a 
4 



o 

II 

> 




f N l° ) 


- W 1 




D 

o 


CD 

II 

o 




Mj'Va 






o 


V 2 




N 2 


f w 2 


! 


Q 2 ° 


1 i 


a/t O / 



v 3 

-W 3 
a0 o 

v 4 

-w 



+ 'M 2 °/a v (6-1) 



= 0 



a0 4 = 0 



N 3 

q 3 ° 

M 3 °/a 

N 4 ° 

q 4 ° 

o , 

M, /a 



where the zero entries in the and matrices are 

obtained from the symmetry characteristics of the load. The 
initial forces in [x°] are derived in Appendix B and presented in 
Table BI. They are transposed to the left-hand side of Equation 
(6-1) to give 
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N 



M 2 /a 

0 

0 

0 

0 

0 

0 

N„ 



M 4 /a 



12709678^ 

22182251 

01944670 

12707678 

22182251 

01944670 



0 \ 



-w 



-w . 



3. 0" 



El r l 

= -3 [k J 12 X 12 < 



^ ) 



(6-2) 



-w 3 

a0 3 

0 

- w A 



The reduced matrix equation is created from Equation ( 6 - 2 ) by 
striking out the rows of the equation and columns of 
12 x 12. The reduced equation with ^Kpj an 8 x 8 matrix 
appears in Table VIII. 

A computer program, "BETA 30, " a copy of which appears 
in Appendix C, solves Equation (6-3) for ^u^j by using the 
subroutine "DSIMQ. " For this example the output of I'DSIMQ" is 

-0.0596 



w ! = 



w T El 



a 3 P 



v 9 _ v 2 E I 
2 “ 



w 



2 = 



W2EI 

a J P 



0 ?= 3 ^ 

2 a 2 P 



0. 0255 



-0. 0287 



- 0 . 0802 
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REDUCED MATRIX EQUATION FOR CASE 2 LOAD, /3 = 30 DEGREES 



CM 
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cO 
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TABLE IX 



FINITE ELEMENT RESULTS FOR CASE 1 LOAD 



♦ 


: 


ft - 7. 5° 


II 

1 — • 
Ln 
O 


ft = 30° 




W 


-0. 0744 


-0. 0744 

1 


-0. 0744 


0° 


N. 

<t> 


0. 000 


0. 000 


0. 000 




M 

° 


0. 318 


0. 318 


0. 318 


i 




w 


-0. 0334 


-0. 0334 


-0. 0334 


LO 

O 

O 




-0. 250 

| 


-0. 250 

1 


-0. 250 
. 


\ 

i 

i 

' ' 


M 

-J> 


0. 068 

| 


0. 068 


S V 

0. 068 ! 

J 


i 

i 


w 


0. 0364 


0. 0364 


0. 0364 


O' 

o 

o 


■ N ? 


-0. 433 


-0. 433 


-0. 433 






-0. 115 


-0. 115 


\ 

-0. 115 




w 


0. 0683 


0. 0683 


-0. 0683 


o 

o 

O' 


N 


-0. 500 


-0. 500 


-0. 500 




*4 


-0. 182 


-0. 182 


-0. 182 
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TABLE X 



FINITE ELEMENT RESULTS FOR CASE 2 LOAD 



♦ . 




0 = 7.5° 


/3 = 15° 


o 

o 

CO 

II 




w 


-0. 0596 


-0. 0596 




-0. 0596 


0° 


N 

<t> 


-0. 128 


-0. 128 




-0. 128 






0. 190 


0. 190 




0. 190 


. 


w 


-0.0287 


-0. 0287 




-0. 0287 


30° 


% 


-0. 239 


-0. 239 




-0. 239. 




M* 


0. 080 


0. 080 




0. 080 




. w 


0. 0298 


0. 0298 




0. 298 


o 

o 


N" 

4> 


-0. 413 


-0. 413 




-0. 413 






-0. 095 


-0. 095 




-0. 095 




w 


0. 0574 


0. 573 




-0. 0573 


o 

o 

o 




-0. 477 


-0. 477 




-0.477 ' 

j 






-0. 159 


-0. 159 




-0. 159 
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TA BEE XI 



FINITE ELEMENT RESULTS FOR CASE 3 LOAD 



♦ 




/3 = 7. 5° 


tl = is 0 


= 30° I 




w 


-0. 0287 


-0. 0287 


1 

-0. 0287 I 


0° 


** 


-0. 239 


-0.239 


-0. 239 




m 6 


0. 080 


0/080 


0. 080 




w 


*• 

-0. 0149 


1 . 

-0. 0149 


-0. 0149 


30° 


N. 


-0. 271 


-0. 271 ’ ; 

i i 


-0. 271 




; M 6 


0.048 ' 


1 i 

0. 048 


0. 048 1 




w 


0. 0143 


o. 0143 


1 

0. 0143 


60 ° 




-0. 358 


. -0. 358 


-0. 358 






-0. 040 


-0. 040 


-0. 040 




w 


0. 0298 


0. 0298 


-0. 0298 


90°. 


N. 


-0. 413 

1 


-6. 413 

| 


-0. 413 


i ! 

• < 




-0. 095 


-0. 095 


-0. 095 
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CHAPTER VII 

COMPARISON AND DISCUSSION OF THE RESULTS 

In order to provide a standard to which the results of the 
Fourier series and finite element methods can be compared, the 
ring problem is solved analytically for N^ and in Appendix A 
for the three load conditions. The results of this analytical 
solution are non -dimensionalized and are entered in Tables XII, 

XIII and XIV with the results from the Fourier series and finite 
element methods. The Fourier series values in these three 
tables are the converged values of w, N^ and while the finite 

element values are the same quantities obtained by using ft = 7. 5 
degrees. A comparison of the results entered in Tables XII, 

XIII and XIV shows essentially zero error. This is well within 
normal engineering tolerances for structural analysis and indicates 
that under proper conditions either of the two methods is capable 
of satisfactory accuracy. 

Since accuracy alone does not provide a sufficient basis for 
choosing one method over the other, the convenience or ease of 
usage of each method must be examined. Ease of usage must be 
considered within the framework of the load applied, the accuracy 
desired, the equipment available, and the assumptions made. A 
comparison of the ease of usage of the two methods within the 
framework of the load applied reveals that the Fourier series 
method has distinct advantages over the finite element method when 
solving the ring problem for a load that is applied over the surface 
of the ring in such a manner that P^ is a function that is easily 
evaluated for changing m and ([). There are two reasons for this. 
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TABLE XII 



RESULTS FOR CASE 1 LOAD 







Fourier Series 


Finite Element 


A nalytical 




w 


-0. 0744 


-0. 0744 


- 




0° 


N 


0. 000 


0. 000 


0. 000 








0. 318 


0. 318 


0. 318 






w 


-0. 0334 


-0. 0334 


- 




30° 




-0. 250 


-0. 250 


-0.250 








0. 068 


0. 068 


0. 068 






w 


0. 0364 


0. 0364 


- 




0' N 

o 

o 




-0. 433 


-0. 433 


-0. 433 








-0. 115 


-0. 115 


-0. 115 






w 


0. 0683 


0. 0683 


- 




o 

o 

O' 




-0. 500 


-0. 500 


< 

-0. 500 








-0. 182 


-0. 182 


-0. 182 
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TA BLE XIII 



RESULTS FOR CASE 2 LOAD 



♦ 




Fourier Series 


Finite Element 


A nalytical 




w 


-0. 0596 


-0. 0596 


- 


0° 


N, 

t 


-0. 128 


-0. 128 


-0. 128 




M. 

<p 


0. 190 


0. 190 


0. 190 




w 


-0. 0287 


-0. 0287 


- 


30° 




-0. 239 


• -0.239 


-0. 239 






0. 080 


0. 080 


0. 080 




w 


0. 0298 


0. 0298 


- 


60° 


S » 


-0. 413 


-0. 413 


-0. 413 






-0. 095 


-0. 095 


-0. 095 




W t ■ 


0. 0574 


0. 0574 


- 


o 

o 

o . 


’ V- 


-0. 477 


-0. 477 


-0. 477 






-0. 159 


-0. 159 


-0. 159 

- - 
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TABLE XIV 



RESULTS FOR CASE 3 LOAD 



♦ 




Fourier Series 


Finite Element 


A nalytical 




\v 


-0. 0287 


-0. 0287 




0° 




-0. 239 


-0. 239 


-0. 239 




s 


0. 080 


0. 080 


‘ ' 0/080 




w 


-0. 0149 


-0.0149 


, 


30° 


N. 

t 


-0. 271 


-0.271 


-0 .271 ' 






0. 048 


0. 048 


0. 048 

i 




w 


0. 0143 


0. 0143 


- 


O' 

o 

o 




-0. 358 


-0. 358 


-0.358 




Mj, 


-0. 040 


-0. 040 


-0. 040 




w 


0. 0298 


0. 0298 




90°’ 


N, 


-0. 413 


-0. 413 


-0. 413 






-0. 095 


-0. 095 


-0. 095 
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The first is, that the series is easily evaluated when r" is a 
simple function. The second is that for distributed loads the 
finite element method required initial fixed end conditions for 
those elements containing the loads. Unless these initial end 
conditions are available from a handbook or other source, they 
must be calculated using energy principles or some other method. 
As can be seen in Appendix B, calculation of these initial end 
conditions may be laborious and time consuming and may re- 
quire a high degree of initial accuracy for the results to be 
useful. Even if the equations for calculating the initial end 
conditions are available in a handbook, they may be lengthy 
and may also require a high degree of initial accuracy. 

A situation in which the finite element method is superior, 
in the framework of the load applied, is one where the load con- 
ditions lead to either |x°^ — 0 initial end conditions or fixed end 
conditions that are available and easily calculated. This is 
especially significant if several such loading conditions are 
to be considered. All that is required for a solution in such 
a situation, once a computer program has been written, is to 
change the program inputs to reflect the changes in[x} , 

|x°] and |uj caused by the load change and to run the program 
for each loading condition. If the Fourier series method is 
used in this situation, each load change requires an integration 
to obtain P , and a different series must be calculated for each 
loading condition. 

As pointed out earlier in this chapter, both methods are 
capable of giving excellent accuracy under certain conditions. 
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For the Fourier series method, accuracy is a:i>_ f^d only b ^ the 
number of terms included in the summation, pro; ess. For dis- 
placements, Tables U, 111 and IV show that very few terms are 
required for a high degree of accuracy especially when a distri- 
buted load is applied. Using only 4 terms in any of the summations 
for w, including that for the point load, gives a maximum error 
of 0. 02 per cent. For and the situation is not quite as 
favorable, but at the worst, 4 terms give a maximum error of 
3. 5 per cent. If this degree of accuracy is satisfactory, the 
Fourier -enes method must be considered easier to use than the 
finite eleinent method. 

For the finite element method, ail element sizes considered 
gave essentially the same results for all three load conditions. 
However, for load conditions o*her than the ones considered here, 
the analyst may be required to use elements with /3 less than 7. 5 
degrees. There appear to be two limiting factors cn the minimum 
size of elements that can be used. The first limiting factor is 
computer capacity. As the eleme nts are made smaller, more 
nodes are created, and hence larger matrices are required. When 
dealing with full scale structures, the limit of computer capacity 
can be reached prior to reaching the optimum element size for 
accuracy. The other limiting factor is related to computer 
capacity but arises from a characteristic of the curved element 
stiffness matrix rather than from the number of elements used. 

The characteristic referred to is the tendency for some of the 
numbers in the stiffness matrices, particularly on the main diago- 
nal , to have a greater disparity ip. size as (Z is decreased. Tables 

V, Viand VII show this to a marked degree. For instance, the 

5 } 



ratio of the stiffness coefficient (6, 6) to coefficient (1, 1) in 
j^K^for /3 — 3 0 degrees is 17/17,243, whereas the ratio of 
the same two elements is 69/ 18, 665, 685 for ft = 7. 5 degrees. 
Since the computer retains only a certain number of significant 
figures while manipulating these matrices, a point will be 
reached where round-off error in computer operations will 
negate any gain in accuracy obtained by using smaller elements. 
This problem is alleviated to a degree by using double-precision 
techniques, as used in this thesis; but again since double- 
precision operations require more computer storage space, the 
computer capacity itself places a limit on this approach. 

The equipment available to the analyst has a significant 
effect on the ease of usage of the two methods. Indeed, in the 
case of the finite element method, lack of a digital computer 
almost prohibits use of the method on the basis of effort required. 
Even if a digital computer is available, a desk calculator and a 
set of mathematical or trigonometric tables with up to ten signifi- 
cant figures is desirable if the nature of the load is such that 
fixed end conditions must be determined. On the other hand, 
while a desk calculator, and even a digital computer, may be used 
to advantage in the Fourier series method, neither are required.'’ 
A slide-rule, a table of integrals, and normal trigonometric 
tables are sufficient equipment for application of the Fourier 
series m ethod. 

Finally, the ease of usage must be considered within the 
framework of the theory employed. The assumptions made when 
applying the two methods of solution to the ring can have an effect 
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on both the accuracy obtained and the e, us ig. . .. 1 

governing differential equations are derived, and the Fourier 

series solution developed through Equation (4-6a), (4-6b) and 

(4-6c) without making any assumptions regarding the extension 

of the middle surface of the ring. If the assumption is made 

during the derivation of these differential equations that the 

middle surface is inextensional, that is £ ^ 0 when z - 0, then 
d v 

^ — : -w and the m r 0 , or extensional, terms do not appea r in 
the expressions for w and M^. * As seen in Table I, no large 
error is introduced in w and by this assumption. However, 
if this inextensional form of w is used in Equation (2-6a) to 
obtain N^, the resulting expression is missing the leading, or 
m = 0, term, and the values for N . calculated from this shortened 
expression are in unacceptable error. This problem may be 
circumvented by using the equilibrium equation, Equation (2-2a), 
to calculate N^. ** The procedure used is to substitute the inex- 
tensional form of M^, obtained from Equation (2-6b), into 
Equation (2-2a) and integrate. The constant of integration is 
obtained from equilibrium considerations at one of the planes of 
symmetry. The resulting expression for is identical to 
Equation (4-4b) and thus gives accurate results. On the other 
hand, in the finite element method, an error will occur in the 



* The governing equations for the modes m > 2 identically 
satisfy the inextensional assumption. Thus, these nodes are 
not affected by the assumption. 

** This is analgous to determining from the equilibrium 
equation, Equation (1-la). 
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solution for w if the inextensional assumption is m^ide and 
distributed loads are applied over almost all the ring surface. 

The reason for this is that a distributed load acting on the surface 
of a curved element has a greater tendency to compress or ex- 
tend the middle surface than does a concentrated load. Errors 
are caused if the energy of this middle surface compression or 
extension is neglected, as it is in the derivation of the curved 
element stiffness matrix. Unlike the Fourier series method, 
there is no means of introducing this extensional behavior into 
the finite element method once the stiffness matrix is derived. 

The only recourse is to select a t/a ratio and re-derive the 
curved element stiffness matrix including the axial strain energy 
term in the calculations. 
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CHAPTER VIII 



CONCLUSIONS AND RECOMMENDATIONS 

The primary advantage of the Fourier series method for 
analyzing ring structures is that it gives reasonably accurate 
answers with a minimum amount of effort using a slide rule 
and normally available integral and trigonometric tables. 

The primary advantage of the finite element method lies 
in its ability to accurately solve several different loading con- 
ditions with a minimum amount of effort. The finite element 
method has the disadvantage of requiring the calculation of 
fixed end conditions for certain types of loads. 

A .recommendation for further work is the derivation of a 
stiffness matrix for a curved element that includes the effects 
of axial strain energy. Also recommended is a matrix error 
analysis on both the extensional and inextensional forms of the 
curved element stiffness matrix for various values of the 
central angle. 
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APPENDIX A 



ANALYTICAL SOLUTION 



The ring problem is solved analytically for N^ and for 
all three load conditions to obtain a standard to which the Fourier 
series and finite element results can be compared. The ana- 
lytical solution is based on the symmetric character of the loads 
applied, the principle of superposition and Castigliano ! s Theorem. 
The use of superposition in the solution of thin shell problems is 
discussed in detail in Reference 5. Figure A 1 (I) shows the 
entire ring with an applied load of P pounds uniformly distri- 
buted over each of the two surface areas subtended by the angle 
2 <=< . Figures A 1 (II), (III) and (IV) show the portions of the 
ring that are superimposed to obtain a solution. 

Considering the quadrant of the ring for 0 — (j) If tf/2, the 
boundary conditions plus equilibrium for Figure A 1 (II) require 



that 




(A - la ) 



and 




(A -lb) 



for 0— § — . Similarly 



( N ^)ni ~ H 2 COS ^ 



(A -2a) 



and 





for Figure A 1 (III). 
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PAbO 




FIGUEE A 1 

EING SEGMENTS USED IN ANALYTICAL SOLUTION 
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In Figure A 1 (IV ) 



Psin cx sin ({) 



^ N <f> W 2 °< 



(A -3a) 



and 



, a PsinoC ( 1 -sin <t)) 

< M ^IV = M l + 2^ “• 

are written for c< < (j) < IT / 2. Since equilibrium gives 



M 



1 



rr - a cos o< 



aP ( 1 -sino< ) 



can be substituted in the expression for (M^py to obtain 
When (|) = 



\ \ a TT _ . aP (l-sinc*sin (D) /A Q1 v 

(M^iy - M 2 " a H 2 cos c< + g - - ■“ (A -3b) 



Mj 



M. 



> ' " -2 
so that Equation (A -3b) yields 



H. 



Pcos° <. 

2 



This expression is substituted into Equations (A -2a) and (A -2b), 
and the resulting equations are superimposed on Equations (A -la) 
and (A -lb) to obtain 



Nj 



(C0So< COS (p - 1 ) 



(A - 4a ) 



and 



M 



* 



^ ^ (cos (J) - cos®< ) (A - 4b) 



for Of (ff 



Substituting the same expression into Equation (A -3b) yields 

a P sin 






m 2+ 



(sin^ - sin (j))- 



(A -5) 



where ^ (p ^ TT / 2 . 

Castigliano's Theorem is applied in the form of 
1 ^ ^ r: 0 to obtain the internal moment The total 

strain energy, U t> of the three segments is written as 






U IIa r U Hb~ t U IIIa ‘ u IIIb+ u IVafUlVb 
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where subscripts a and b stand for axial and bending strain 
energy respectively. The terms Ujj a ,’ Uj^, an ^ Ujy a neec * 

not be considered when applying Castigliano's Theorem to obtain 
since these terms are either zero or are not functions of M 2 . 

>>Ut 



Therefore 



0 bec 



om es 



V_ 



+■ 

Z El 



-1T/z x 



Ztl 



= 0 



(A -6) 



where the expressions for are obtained from Equations (A -4b) 
for the first integral and (A -5) for the second integral. Performing 
the partial differentiation and integration with the appropriate 
limits as indicated in Equation (A -6) gives 




aP 

IT 



f r my\ <=< 



2. Ot 

9 




(A “7 ) 



Using this expression for M in Equations (A -4b) and (A -5), the 

Cd 

complete set of equations for and is written. 



N 



<t> 



P 

2 *< 



cow cos -l 
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1 



1 



CO^c*. cos dp 

?-<=< 



M* = aP 



IT 



2=* 



for 0 ^ <{> r °< , and 




P S \ v\ <=*- s t n <J> 




s m c< s m (t> 
Z°( 



foro< ^ (j) ^ TT / 2 . , These equations are non-dimension- 

alized and evaluated for (J> = 0, 30, 60 and 90 degrees for all 

three case loads. The results are entered in Table A I. 
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TABLE A I 



ANALYTICAL SOLUTION RESULTS 



<p 




Case 1 Load 


Case 2 Load 


Case 3 Load 




N. 

<t> 


0. 0 


-0. 128 ' 


-0. 238 


0° 




0. 318 


0. 190 


0. 080 




N* 


-0. 250 


-0. 239 


-0. 270 


o'*. 

, o 
* 












*♦ 


-0. 068 


0. 080 


0. 048 




N<t> 


-0. 433 


-0. 413 


-0. 358 


o 

o 












si 

-e- 


-0. 115 


-0. 095 


-0. 040 






-0. 500 


-0. 477 


-0. 413 


o 

o 

o 














-0. 182 


-0. 159 


-0.095 
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APPENDIX B 



SOLUTION FOR INITIAL FORCES 



The techniques used to obtain the initial fixed end moments 

and forces required for the finite element solutions to Case 2 and 

Case 3 loads involve superposition and Cas tigliano’s Theorem. 

Figure B 1 shows the finite element, the three loading conditions, 

the boundary forces and moments, and the notation used. The 

superposition of the forces and moments in Figure B 1 (II) and 

B 1 (III) is equivalent to the forces and moments in Figure B 1 (I). 

The sum of the rotations of the ends of the element in Figure 

B 1 (II) and (III) must be zero. Thus, 0jj + = 0. 

However, since 0 is zero under the end conditions and loading 
II 

of Figure B 1 (II), 0 jjj = 0. The analogous requirement for the 
horizontal displacements of the ends of the element of Figure 
B 1 (II) and (III) leads to 



from an expression for the strain energy of the element in Figure 




(B-l) 



where is the horizontal component of the known radial dis- 

placement of (II), and Sjjj is the horizontal displacement obtained 



B 1 (III) using Castigliano 1 s Theorem. 



direction towards the centerline, and 




positive in the 



positive in the 



direction away from the centerline. 
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2<* a b ? 





CQ 

W 

« 

P 

a 

i— i 



RING SEGMENTS USED IN FIXED END SOLUTION 



Letting cs be the half angle of the element and <|) be the 
angle at any point along the element, the moment at ({) is 

M $ = - aHrrt 1 cos 0 - c o s®<] 

And the axial force at Q is 



= H m cos 

The strain energy, from both axial extension and bending 
of the total, element is 



U t = 2- 



r.'V 



Z El 



r' 

l 

i 2 AE 






(B-2) 



Using Equation (B-2) in Castigiiano 1 s- Theorem gives 



2 { - ^ ^ (B 

w *- 

Substituting Equation (B-2) into Equation (B-3), differentiating, 
integrating and applying limits yield 



C 

^nr - 



EX 1 



a 1 - -Ljs\v\o<-<* 



3 vn 

h- 






2 . 

cos ©< 



Q p, nr 
AE. 



o< 

X 



(B-4) 



SW 2 
H- 



J 
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The horizontal displacement 



is 



II 




a P s \Y\<=< 
Z«* V>-ft £ 



(B-5) 



where f is given by 



1 

CASE LOAD 


2 


3 


(3= 7.5° 


Mi 

II 

00 


f = 16 


P = 15. 0° 


II 


II 

00 


J = 30. 0° 


l-h 

II 

ro 


f = 4 

> 



The factor f is required in order to give the proper portion of the 
total load P over each of the different segment sizes. Substi- 
tuting Equations (B-4) and (B-5) into (B-l) leads to 



— 

BI 111 



— — (" s \ °< — Os cos <=^*1 



3 sm 
H 



+ c< COS c< 



+ AE H nr 



C* MVI t<=< 

— + 

2 - H- 



a P s 



ur\«* 






The above equation can be simplified to. give 



a T~ 

12 . It) h 



nr 



i_ 



2 

^ S \Y\ o<-ccCOS°s] + ~ 



3 S \y\Z<* 



+ ex' c o S =< 



+- H 



nr 



<x sm2°<- 

— -l- 

X h 



P S t r\ 



2.<=< f 



(B-6) 



Thus, Hj.jj.can be determined from Equation (B-6) for a given case 
load and a specified central angle for the element, /3 , once t/a is 
chosen. The requirement to select a specific t/a ratio in order to 



68 



solve Equation (B-6) means that the axial strain energy is in- 
cluded in the expression. Failure to include the axial strain 
energy by making the assumption made in Reference 2 that it can 
be neglected for curved elements with small t/a ratios leads to 
gross errors. This is because the assumption made in Refer- 
ence 2 is not valid for curved elements with distributed loads 
over large portions of the element. 

The equation for the total normal forces N° acting on the 
ends of the element shown in Figure B 1 (I) is 



Nnt- H m cosc< 



or 




(B-7) 



since 



-f 



The component of in the radial direction is the shear 

force, Q°, at the ends of the ring element. Thus 



Q° = Sffl« 



(B-8) 



In order to obtain M°, Castigliano 1 s Theorem, in con- 
junction with the requirement 0 T j T = 0, is used to write 
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Differentiation and integration of this equation give 




$>\V\oc - 



COS o< 



(B-9) 



Experience with Equations (B-6) through (B-9) shows that the 
axial strain energy must be included if accurate answers are de- 
sired for N°, Q° and M°. Thus, t/a = 1/10 is chosen, and 
Equations (B-6) through (B-9) are used to prepare the entires 
in Table B I for N°, and M° for Case 2 and 3 loads and /3 = 30, 

15 and 7. 5 degrees. In evaluating these equations, trigonometric 
functions accurate to 5 decimal places do not give accurate 
answers for N°,*Q° and M° for the values of used in this 

thesis. Consequently the entries in Table B I are obtained by 

\ 

using trigonometric functions of 10 decimal place accuracy. 
Solving equations such as (B-6) using 10 decimal places is ob- 
viously beyond the capacity of either slide rule or practical hand 
calculations. Even with an electronic desk calculator, solving 
such equations is laborious and time consuming. 



70 



TABLE B I 



FIXED END RESULTS 



CASE 2 LOAD 





7. 5° 


15° 


3 0° 


N° 


0. 00181929P 


0. 01286764P 


0. 12707678P 

i 


Q° 


0. 06247 01 5P 


0. 12402482P 


0. 2218225 IP 


M° 


0. 00136328Pa 


0. 00541 780Pa 


0. 0 19446 70Pa 



CASE 3 LOAD 





7. 5° 


15° 


30° 


N° 


0. 00090964P 


0. 00643382P 


0. 06353839P 


Q° 


*0. 03 123508P 


0. 06201237P 


0. 1 1091 125P 


M° 


0. 00068 164Pa 


0. 0027089 OPa 


0. 00972335Pa 
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A PPENDIX C 



FORTRAN IV PROGRAMS 
TABLE C I 

DESCRIPTION OF SUBROUTINE "DSIMQ" 



PURPOSE 

OBTAIN SOLUTION OF A SET OF SIMULTANEOUS LINEAR 
EQUATIONS, AX=B 

USA GE 

CA LL "DSIMQ" (A , B, N, KS) 

DESCRIPTION OF PARAMETERS 
A AND B MUST BE R EA L*8 

A - MATRIX OF COEFFICIENTS STORED COLUMNWISE. 
THESE ARE DESTROYED IN THE COMPUTATION. 

THE SIZE OF MATRIX A IS N BY N. 

B.- VECTORS OF ORIGINAL CONSTANTS (LENGTH N). 

THESE ARE REPLACED BY FINAL SOLUTION VALUES 
VECTOR X. 

N - NUMBER OF EQUATIONS AND VARIABLES 
KS -OUTPUT DIGIT 

0 FOR A NORMAL SOLUTION 

1 FOR A SINGULAR SET OF EQUATIONS 

R EMA R KS 

MATRIX A MUST BE GENERAL. 

IF MATRIX IS SINGULAR, SOLUTION VALUES ARE MEAN- 
INGLESS. AN ALTERNATIVE SOLUTION MAY BE OB- 
TAINED BY USING MATRIX INVERSION (MINV) AND 
MATRIX PRODUCT (GMPRD). 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
NONE 

METHOD 

METHOD OF SOLUTION IS BY ELIMINATION USING 
LARGEST PIVOTAL DIVISOR. EACH STAGE OF ELIMINA- 
TION CONSISTS OF INTERCHANGING ROWS WHEN NECES- 
SARY TO AVOID DIVISION BY ZERO OR SMALL ELEMENTS 
THE FORWARD SOLUTION TO OBTAIN VARIABLE N IS 
DONE IN N STAGES. THE BACK SOLUTION FOR THE 
OTHER VARIABLES IS CALCULATED BY SUCCESSIVE 
SUBSTITUTIONS. FINAL SOLUTION VALUES ARE 
DEVELOPED IN VECTOR B, WITH VARIABLE 1 IN B (1), 

VARIABLE 2. IN B(2) VARIABLE N IN B(N). 

IF NO PIVOT CAN BE FOUND EXCEEDING A TOLERANCE 
OF 0. 0, THE MATRIX IS CONSIDERED SINGULAR AND 
KS IS SET TO 1. THIS TOLERANCE CAN BE MODIFIED 
BY REPLACING THE FIRST STATEMENT. 
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TABLE C II 



LIST OF SYMBOLS USED IN SUBROUTINE "EELM" 



Computer 
Coded Name 

DBETA 

DESM 

DESM (I, J) 

DLA 

DLB 

DLC 

DLD 

DLE 

DBA 

DBB 

DBC 

DBD 

DBE 

DBF 

DBG 



Definition 

/3 , element central angle. 
Element stiffness matrix . 

Stiffness influence coefficient in 

[H • 

a (corresponding notation in 
b Reference 3). 

c 

d 



e 

A 

B 

C 

D 

E 

F 

G 
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TABLE CIII 



SUBROUTINE "RELM" 



SUBROUTINE REL M ( OB ETA , DESM f OL A , DIB, DLC. , OLD , DL E , OB A , 
1DBB, DBC.DBD.DBE, DBF. DBG) 

OOUBLE PRECISION DLA, DLB , DLC , OLD, DL E, DBA , DB8 , DBC , DBD, 
1DBE,DBF, DBG 
DOUBLE PRECISION DBET A 
REAL*8 DESMI 6,6 ) 

DLA=DBETA-D$IN(DBETA) 

DL B=DCOS I DBETA ) ■*■ I ( DS INI DBET A) ) **2 ) / 2 . ODO- 1 . ODO 
DLC= ( 3 .ODO*DBETA ) /2 .0 00-2 .ODO*DS I N( DBETA ) ♦ 
1(DSIN(2.0D0*DBETA) 1/4.0 DO 
DLC=DBETA/2.0D0-(DSIN(2.0D0*DBETA) J/4.0D0 
DLE=DCOS(DBETA)— l.ODO 
DBA= I DLE **2 ) /DBETA-DL D 
DBB=DLB- ( DL A*DLE )/ DBETA 
DBC=DLA*DLD-DLB*DLE 
DBD=(DLA**2) /DBETA-DLC 
DB E=DLC*DL E-DL A*DL B 
D8F=(DLB**2)-( DLC*DLD ) 

DBG=DLB*(DLB-(2.0D0*DLA*DLE)/DBETA)+DLC* 

1I-DL0+I (DLE**2) /DBETA) )+( ( DLA**2 )*DLD) /DBETA 
DE$M(1,1)=DBA/DBG 
DESM(1,2)=DBB/DBG 
DESM(l,3)=DBC/( DBETA* DBG) 

DESM(2,1)=DBB/DBG 

DESM(2,2)=DBD/DBG 

DESM(2,3)=DBE/(DBETA*DBG) 

DESM(3,1 )=DBC/(DBETA*DBG) 

DE$M(3,2)=DBE/(DBETA*DBG) 

DESMI 3 ,3 )=DBF/ I DBET A* DBG) 

DESM 1 4,1 )=-( DBA*DCOS I DBETA )+DBB*DS INI DBETA) ) /DBG 
DESM 1 4,2)=- I DBB*DCOS I DBET A ) +DBD*DS INI DBETA) ) /DBG 
“ I DBC*DCOS( DBET A)+DBE*DSIN( DBETA) )/ 



DESM(4,3 )=- 
l I DBETA*DBG ) 
DESM(5,1)= 
DESM(5,2)= 
DESMI 5 ,3 )= 

1 I DBE T A*DBG ) 
DESMI6, l )= 



( DBA* DS I N I DBET A ) -DBB*DCOS I DBET A) ) /DBG 
(DBB*D$IN(DBETA )-DBD*DCOS I DBETA )) /DBG 
I DBC*DS INI DBET A J -DBE*DCOS I DBET A ) )/ 



(DBA*! DCOSIDBETAJ-l.ODO) + DBB* 
1 1 DS I N( DBET A) ) -DBC/ DBETA ) /DBG 
DESM ( 6,2)= (DBB*(DCOS( DBET A)-l . ODO > +DBD* 
IDS INI DBETA)- DBE/ DBETA) /DBG 
DESM I 6,3)= I DBC* I DCOS I DBET A ) -1 . ODO ) +DBE* 
IDSINI DBETA)- DBF)/! DBETA*DBG) 

DESM(4,4)=DESM( 1.1 ) 

DESM 1 4, 5 )=-DESM( 1,2) 

DESM(4,6)=DESM( 1 ,3) 

DESM(5,4)=-0ESM(2,1) 

DESM(5,5)=DESM(2,2) 

DESM(5,6)=-DESM(2,3) 

DESM(6,4)=DESM(3,1 ) 

DESM(6,5)=-DESM(3,2) 

DESM(6,6)=DESM(3,3) 

DESM(1,4)=DESM(4,1 ) 

DESM(1,5 )=DESM(5,1 ) 

DESM(1,6)=DESM(6,1 ) 

DESM(2,4)=DESM(4,2) 

DESM(2,5)=DESM(5,2) 

DESM(2,6)=DESM(6,2) 

DESM(3,4)=DESM(4,3) 

DESM(3,5)=DESM( 5,3) 

DESM(3,6)=DESM(6,3) 

RETURN 

END 
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TABLE C IV 

LIST OF SYMBOLS USED IN 
PROGRAM "BETA 30" 


Computer 
Coded Name 


Definition 


A 


r i 

Reduced stiffness matrix K p . 


A (I,J) 


Stiffness influence coefficient 
in K^>j. 


B (I) 


Prior to calling "DSIMQ" B(I) = 

f {Xl - fx°l} . 


B (I) 


After calling "DSIMQ" B (I) - [up] 


DBETA 


(i , element central angle. 


DESM 


Element stiffness matrix [k^j . 


DESM (I, J) 


Stiffness influence coefficient 
in r K^j . 


PI 


IT 


DM (I) 


Moment at node (I). 


DN (I) 


Normal force at node (I). 
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TABLE CZ 



PROGRAM "BETA30" 



PROGRAM BETA 30 

DOUBLE PRECISION DBET A. P I , DLA , DL B, DLC , OLD, DLE , DB A , 
1DBB»DBC,DBD,DBE»DBF,DBG 

91 FORMAT C 1H0 , 1 P6D20 . 10 ) 

92 FORMAT ( 1H0 , 1P5D19 . 10 ) 

93 FORMAT ( 1H0 , 1P7D18 . 10 ) 

94 FORMAT < 1H0, 1PD20. 10 ) 

11 FORMAT ( 1H0, 1P5D20.10) 

12 FORMAT ( 1H0,1PD20. 10) 

13 FORMAT ( 1H0, 1P6D20.10 ) 

14 FORMAT ( 1H0, 1P9D15.7) 

15 FORMAT ( IH1 , 'DISPLACEMENTS* ) 

16 FORMAT (1H0, 'FORCES AND MOMENTS') 

17 FORMAT ( 1H0» 'CASE 2 BETA=30 DEGREES') 

22 FORMAT C 1H0, 1PD20.10) 

23 FORMAT (12) 

24 FORMAT ( 1H0 , 1 PD 20 . 1 0 ) 

PI =3 .141 592653589793 DO 
DBET A=P I /6 . ODO 
WRITE (6,94) DBETA 
RE AL*8 DESM ( 6,6 ) 

CALL RELM ( DBETA, DESM, DLA, DLB, DLC, DLD, DLE, DBA, 
1DBB.DBC, DBD,DBE,DBF,DBG) 

WRITE (6,92) DLA, DLB, DLC, OLD. DLE 
WRITE (6,93)DBA,DBB,DBC#DBD,DBE,DBF, DBG 
WRITE (6.91) ((DESM(I,J),J=1,6),I=1,6) 
F0RBETA=30DEGREES 
RE AL*8A (8.8) 

DO 110 1=1,8 
DO 111 J=1 , 8 
A( I , J)=O.ODO 
111 CONTINUE 
110 CONTINUE 

A ( 1, 1)=DESM( 2,2) 

A ( 2 , 1 ) =DESM( 4,2) 

A( 3,1) =DESM( 5,2) 

A(4,1)=DESM(6,2) 

A( 1,2) =DESM( 2,4) 

A( 2 , 2 ) =DESM ( 4,4)+DESM(l,l) 

A( 3,2)=DESM(5,4)+DESM(2,1) 
A(4,2)=DESM(6,4)+DESM(3,1) 

A ( 5,2) =DESM( 4,1) 

A(6,2)=DESM(5,1) 

A( 7 , 2 ) =DESM ( 6,1) 

A( 1, 3)=DESM( 2,5) 

A( 2,3)=DESM(4,5)+DESM( 1,2) 
A(3,3)=DESM(5,5)+DESM(2,2) 

A( 4,3)=DESM(6,5)+DESM(3,2) 

A( 5,3)=DESM(4,2) 

A( 6,3) =DESM( 5,2) 

A( 7 , 3) =DESM (6,2) 

A( i,4)=DESM( 2,6) 

A(2,4)=DESM(4,6)+DESM(1,3) 

A(3,4)=DESM(5,6)+DESM(2,3) 

A(4,4)=DESM(6,6)+DESM(3,3) 

A( 5, 4) =DESM ( 4,3) 

A(6,4)=DESM(5,3) 

A(7,4)=DESM(6,3) 

A( 2,5)=DESM(1,4) 

A(3,5)=DESM(2,4) 

A(4,5)=DESM(3,4) 

A( 5,5)=DESM(4,4)+DESM( 1,1) 

A(6 f 5)=DESM(5,4)+DESM(2,1) 
A(7,5)=DESM(6,4)+DESM(3,1) 

A( 8» 5 )=DESM( 5,1) 
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A( 2,6)=DESM( 1,5) 

A ( 3,6)=DESM(2,5) 

A(4,6)=0ESM(3»5) 

A ( 5,6)=DESM( A, 5 )+DESM< 1,21 
A( 6,6)=DESM( 5,5 )+0ESM(2,2) 

A( 7,6) =DESM( 6,5) +DESM(3,2) 

A ( 8,6)=DESM< 5,2) 

A( 2,7)=DESM( 1,6) 

A ( 3 , 7 ) =DESM ( 2,6) 

A ( 4 , 7 ) =OESM( 3,6) 

A< 5,7)=DESM(4,6 )+DESM(l,3) 

A( 6,7)=DESM( 5,6 )+DESM(2,3) 

A( 7,7)=0ESM(6,6)+0ESM(3,3) 

A ( 8 » 7 ) =DESM ( 5,3) 

A( 5,8)=0ESM( 1,5) 

A( 6,8)=DESM( 2,5) 

A(7,8)=DESM(3,5) 

A( 8,8)=DESM( 5,5) 

21 FORMAT( 1H0 , 1P8D16. 8) 

WRITE (6.21) ( ( A( I, J), J=l,8) ,1=1,8) 

REAL *4 R ( 8 * 8 ) 

00 301 1=1,8 
00 302 J=1 , 8 
R( I , J)=A( I , J) 

302 CONTINUE 
301 CONTINUE 
REAL*8B( 8) 

00 112 1=1,8 
8 ( I ) =0.000 
112 CONTINUE 

C B LOAD MATRIX FOR CASE 2. BET A=30 OEG . 

B( 1)=. 2218225100 
8( 2) =. 1270767800 
B(3)=. 2218225100 
B( 4)=-. 0194467000 
WR ITE ( 6 , 22 ) ( B( I ) , 1 = 1 ,8 ) 

N=8 

CALL DSIMQ ( A, B. N, KS ) 

WRITE(6,23) (KS) 

WRITE (6,12) (B( I ) ,1=1,8) 

WRITE (6,15) 

WRITE (6,12) B(l) 

WRITE (6,12) B( 3 ) 

WRITE (6,12) B ( 6 ) 

WRITE (6,12) B ( 8 ) 

DOUBLE PRECISION ONI , 00 1 , DM 1 , DN2 , 002 , DM2 , 
1DN3,DQ3,DM3, 0N4 , DQ4, DM4 

C FORCES AND MOMENTS FOR CASE 2. BET A=30 OEGREES 

0N1=DESM (1,2)* B(1 )+DESM( 1,4)*B( 2)+0ESM( l, 5)* B(3) 
1+0ESM( 1 ,6) *B(4)+. 1270767800 
DM 1 = DE SM (3,2)* B(l) +0E SM ( 3 , 4 ) *B ( 2 ) +DE SM ( 3 , 5 ) * 

1 B ( 3 ) +DES M ( 3, 6 )*B( 4)-. 0194467000 
DN2=DESM (4,2)* B(l) + OESM ( 4 , 4 ) *B ( 2 ) +DESM ( 4, 5 ) * 

1B( 3)+DESM(4,6) *B(4)-. 1270767800 
CM2= OE SM ( 6 , 2 ) * B ( 1 I +OES M ( 6 , 4 ) * B ( 2 ) ♦DES M ( 6 , 5 ) * 
1B(3)+DESM(6.6)*B(4)+.01 94467000 
0N3=0ESM(1 ,1 )*B( 5)+DESM( 1,2)* B( 6 ) +OESM( 1, 3 )*B{ 7 ) 
1+0ESM( 1,5) *B ( 8 ) 

DM3=DESM(3, 1 )*B( 5 ) +OESM ( 3 , 2 ) * B(6) +DESM(3,3) 
1*B(7)+0ESM(3,5)*B( 8) 

0N4=0ESM(4,1»*B( 5)*DESM(4,2)* B(6) +0ESM(4,3) 
1*B(7)+DESM(4.5)*B( 8) 

0M4=DESM(6,1)*B( 5 ) ♦DESM < 6 , 2 I *B ( 6 ) +DESM ( 6, 3 ) *B ( 7 ) 
1*DESM(6,5)*B(8) 

WRITE (6,17) 

WRITE (6,16) 

WRITE (6,12) DN1,DMI,DN2,DM2,DN3,DM3,DN4,0M4 

STOP 

END 
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